/*
 *	geom_3d_v4b: add more record points - now have up to 10 places where the sound is 
 *		recorded
 *
 *	geom_3d_v4 (a next version of geom_3d_new_v3.c)
 *
 *	geom_3d_v4 -- add feature to output coordinates of marked point in rec_geom file
 *		these points are indicted by "M"
 *
 *	this is version 4 - for now just add a small change to write data on the
 *	solid part of the instrument to a file geometry_instrument.dat
 *	and not the boundary points (as in geometry_solid.dat which is still created)
 *
 *	this is version 3 -- allows for more than 3 layers (which was the limit in
 *		previous verisons) and also gives more control over the flow speed
 *		in the initialization region
 *	
 *	this version will enable me to deal with a transverse flute
 *
 *	set up a nonuniform grid in 3d for use in recorder8 and later programs
 *
 *	usage: geom_3d_v3 <in_file> <recorder_geometry.dat> <output_directory> 
 *
 *	where in_file contains data used to define the nonuniform grid
 *
 *	rec_geom.dat defines the instrument and is similar to (but with a few
 *	more options as) the format used by recorder_v8 and grid_3d_new_v2
 *
 *	output_directory contains all the output files including the vtk files that contain
 *		the overall geometry
 *
 *	the output files are 
 *		x_grid.dat -- contains grid information for the x direction in the form
 *				0	x(0)	dx(0)
 *				1	x(1)	dx(1)
 *				etc
 *				where x(i) is the location along x of grid point i
 *				and dx(i) is the spacing between grid points i and i+1
 *		y_grid.dat -- same info for the grid along y
 *		z_grid.dat -- same info for the grid along z
 *		geometry_solid.dat -- locations of all solid points (recorder and boundary)
 *			in terms of i,j grid locations
 *		geometry_initialize.dat -- locations of the initialization region
 *			in terms of i,j grid locations
 *		geometry_sound.dat -- locations of the points where the sound is to be recorded
 *			in terms of i,j grid locations
 *
 *		and new in v4
 *		geometry_instrument.dat -- locations of all solid points (recorder only)
 *			in terms of i,j grid locations
 *
 *		place recorder_geometry.dat into the output_directory for further ref
 *		
 *		produce eps files showing the precise geometry of the simulation - one set of files
 *		low_res_geom_k.eps is a low res plot and high_res_geom_k.eps shows the detail in the region
 *		where the grid spacing is smallest -- the index k is the layer number as used
 *		in recorder_geometry.dat
 *
 *		the file geometry.vtk is generated for viewing with paraview
 *		and file bare_geometry.vtk is generated for viewing with paraview
 *			(this file assumes a uniform grid)
 *
 *	-------------------------
 *
 *	the way the grid spacing varies for a given direction (i.e., the x direction) is 
 *	defined by several quantities
 *
 *	x_1 <	x_2 <	x_3 <	x_4
 *
 *	x_1: the lowest (starting) index value and dx is large at this end
 *	x_2: dx falls to its minimum value dx_min here and is constant until
 *	x_3: dx is constant between x_2 and x_3. dx grows again for larger values of x
 *	x_4: the largest (final) index value and dx is large at this end
 *
 *	the corresponding array indices are i_1, i_2, i_3, and i_4
 *
 *	add similar code for grid along y direction
 *
 *	this program reads in values for x_1, x_2, .... y_1, y_2 ... from <in_file>
 *
 *	format for in_file is
 *	x_1	x_2	x_3	x_4
 *	dx_min	dx_max	ax
 *	y_1	y_2	y_3	y_4
 *	dy_min	dy_max	ay
 *	z_1	z_2	z_3	z_4
 *	dz_min	dz_max	az
 *
 *	ax, ay, and az are factors that controls how fast dx, dy, and dz vary (higher values give a faster
 *		variation; try a = 20 first)
 *
 *	as noted above, files x_grid.dat and y_grid.dat contain 
 *		i1	x1	dx1	
 *		i2	x2	dx2	
 *		etc
 *	where i1 is the grid index
 *
 *	The files geometry_solid.dat, geometry_initialize.dat, and geometry_record.dat contain
 *	the grid coordinates i,j,k of the solid regions (boundary and recorder), the region where u
 *	is initialized, and the points where the sound is to be recorded
 *
 *	geometry_initialize.dat also contains the scaled values of the fluid speed at each point
 *		on a scale where max = 1.0
 *	
 *	The files x_grid.dat, y_grid.dat, and z_grid.dat  will be used by recorder8 to set up the computational
 *	grid, and the files geometry_solid.dat, geometry_initialize.dat, and geometry_record.dat
 *	will be used to initialize the geometry of the problem
 *
 *	also put in_file into the output directory along with some other files that will not be used
 *		but are useful for record keeping
 *
 *	the file u_init.dat contains info for initializing u in the channel
 */


#include <math.h>
#include <stdio.h>

#include "grid.h"

double exp();
double stretch_factor();

double y_grid[N_Y_MAX];
double dy[N_Y_MAX];
double x_grid[N_X_MAX];
double dx[N_X_MAX];
double z_grid[N_Z_MAX];
double dz[N_Z_MAX];
double ax1,ax2,dx_min,dx_max;
double ay1,ay2,dy_min,dy_max;
double az1,az2,dz_min,dz_max;
double x1,x2,x3,x4;
double y2,y3,y4,y_1;
double z1,z2,z3,z4;
FILE *fp_grid_input;
FILE *fp_x_grid_dat,*fp_y_grid_dat,*fp_z_grid_dat;
FILE *fp_recorder;
double y_max,x_max,z_max;
int i_min,i_max;
int j_min,j_max;
int k_min,k_max;
int i_max_bare,j_max_bare,k_max_bare;
int i_init_min,i_init_max;
int j_init_min,j_init_max;
int k_init_min,k_init_max;
int i_init_min_bare,i_init_max_bare;
int j_init_min_bare,j_init_max_bare;
int k_init_min_bare,k_init_max_bare;

// this needs work
double init_val_bare[N_Y_MAX];
double init_val[N_Y_MAX];
int i_init_bare;
int j_init_bare_start,j_init_bare_end;

// int geometry[N_X_MAX][N_Y_MAX][N_Z_MAX]; /* describes the domain of the calc on NONUNIFORM grid	*/
char geometry[N_X_MAX][N_Y_MAX][N_Z_MAX]; /* describes the domain of the calc on NONUNIFORM grid	*/
                                /* the geometry (including the shape of the recorder)   */
                                /* is specified in recorder_geometry.dat                */
                                /* most geometry-related parameters are set there       */
                                /* OPEN = open space away from walls    */
                                /* SOLID = body of recorder             */
				// LOWER_X_BOUNDARY        boundaies of simulation region
				// UPPER_X_BOUNDARY        
				// LOWER_Y_BOUNDARY       
				// UPPER_Y_BOUNDARY      
				// INIT_POINT              points where u is set

// int bare_geometry[N_X_MAX][N_Y_MAX][N_Z_MAX]; 
char bare_geometry[N_X_MAX][N_Y_MAX][N_Z_MAX]; 
			/* describes the domain of the calc on UNIFORM grid	*/

/*
int i_record_1,i_record_2,i_record_3;
int j_record_1,j_record_2,j_record_3;
int k_record_1,k_record_2,k_record_3;
int i_record_1_bare,i_record_2_bare,i_record_3_bare;
int j_record_1_bare,j_record_2_bare,j_record_3_bare;
int k_record_1_bare,k_record_2_bare,k_record_3_bare;
*/

int i_record[11],j_record[11],k_record[11];
int i_record_bare[11],j_record_bare[11],k_record_bare[11];
int n_record;	// number of points where sound is to be recorded

int i_recorder_corner,j_recorder_corner,k_recorder_corner;
int i_recorder_bare,j_recorder_bare,k_recorder_bare;// bare coordinates of the lower left corner of recorder
		// WARNING - in previous versions (with uniform grids) these were the coordinates
 		// of the lower left corner
int i_init,j_init_start,j_init_end,k_init_start,k_init_end;
double x_recorder_bare,y_recorder_bare,z_recorder_bare;

double dx_bare;	// original "bare" dx for defining the recorder geometry
FILE *fp_log;

char output_directory[200],recorder_geometry[200],input_file[200];
char vtk_directory[200];

FILE *fp_diag;

double v_scale_factor_bare[N_Z_MAX];	/* scaling factors for v_x in each layer that contains an init point in bare coordinates */
double v_scale_factor[N_Z_MAX];		/* scaling factors for v_x in each layer that contains an init point in final coordinates*/

main(int argc,char *argv[])
{
	int i,j;
	char mes[200];

	fp_diag = fopen("diag_info","w");

	if(argc < 4) {
		fprintf(stderr,"usage: geometry_setup <in_file> <recorder_geometry.dat> <output_directory>\n");
		fprintf(stderr,"\tformat for input file is\n");
		fprintf(stderr,"\tx_1\tx_2\tx_3\tx_4\n");
		fprintf(stderr,"\tdx_min\tdx_max\tax\n");
		fprintf(stderr,"\ty_1\ty_2\ty_3\ty_4\n");
		fprintf(stderr,"\tdy_min\tdy_max\tay\n");
		fprintf(stderr,"\tz_1\tz_2\tz_3\tz_4\n");
		fprintf(stderr,"\tdz_min\tdz_max\taz\n");
		fprintf(stderr,"\n");
 		fprintf(stderr,"\ta(x,y,z) (= factor that controls how fast delta varies (try a = 20 first)\n");
 		fprintf(stderr,"\\ttlarger values of ax give a faster variation of dx, etc.\n");
		exit(0);
	}

	fp_log = fopen("log_file","w");
	strcpy(recorder_geometry,argv[2]);
	strcpy(output_directory,argv[3]);
	mkdir(output_directory,0777);
	sprintf(mes,"cp %s %s/",argv[1],output_directory);
	system(mes);
	sprintf(mes,"cp %s %s/",recorder_geometry,output_directory);
	system(mes);
	strcpy(vtk_directory,output_directory);
	

fprintf(stderr,"setting up grid\n");
	set_up_grid(argv[1]);
fprintf(stderr,"reading in and placing recorder\n");
	place_recorder();
	save_geometry();
//	display_geometry(1);
//	make_bare_vtk_file();
	make_vtk_file1();
//	make_vtk_file2();

	check_init_region();

	fclose(fp_diag);
}

set_up_grid(char *in_file)
{
	int i,j,k;
	char mes[200];
	FILE *fp_dx,*fp_dx_display,*fp_x_grid_dat;
	FILE *fp_dy,*fp_dy_display,*fp_y_grid_dat;
	FILE *fp_dz,*fp_dz_display,*fp_z_grid_dat;

	fp_grid_input = fopen(in_file,"r");
	fscanf(fp_grid_input,"%lf %lf %lf %lf",&x1,&x2,&x3,&x4);
	fscanf(fp_grid_input,"%lf %lf %lf",&dx_min,&dx_max,&ax1);
	fscanf(fp_grid_input,"%lf %lf %lf %lf",&y_1,&y2,&y3,&y4);
	fscanf(fp_grid_input,"%lf %lf %lf",&dy_min,&dy_max,&ay1);
	fscanf(fp_grid_input,"%lf %lf %lf %lf",&z1,&z2,&z3,&z4);
	fscanf(fp_grid_input,"%lf %lf %lf",&dz_min,&dz_max,&az1);
	fclose(fp_grid_input);
	
	fprintf(stderr,"%g\t%g\t%g\t%g\n",x1,x2,x3,x4);
	fprintf(stderr,"%g\t%g\t%g\n",dx_min,dx_max,ax1);
	fprintf(stderr,"%g\t%g\t%g\t%g\n",y_1,y2,y3,y4);
	fprintf(stderr,"%g\t%g\t%g\n",dy_min,dy_max,ay1);
	fprintf(stderr,"%g\t%g\t%g\t%g\n",z1,z2,z3,z4);
	fprintf(stderr,"%g\t%g\t%g\n",dz_min,dz_max,az1);
	
	sprintf(mes,"%s/x_grid_spacing",output_directory);
	fp_dx = fopen(mes,"w");
	sprintf(mes,"%s/y_grid_spacing",output_directory);
	fp_dy = fopen(mes,"w");
	sprintf(mes,"%s/z_grid_spacing",output_directory);
	fp_dz = fopen(mes,"w");
	sprintf(mes,"%s/x_grid_display",output_directory);
	fp_dx_display = fopen(mes,"w");
	sprintf(mes,"%s/y_grid_display",output_directory);
	fp_dy_display = fopen(mes,"w");
	sprintf(mes,"%s/z_grid_display",output_directory);
	fp_dz_display = fopen(mes,"w");
	sprintf(mes,"%s/x_grid.dat",output_directory);
	fp_x_grid_dat = fopen(mes,"w");
	sprintf(mes,"%s/y_grid.dat",output_directory);
	fp_y_grid_dat = fopen(mes,"w");
	sprintf(mes,"%s/z_grid.dat",output_directory);
	fp_z_grid_dat = fopen(mes,"w");

	ax2 = (dx_max / dx_min) - 1.0;
	ay2 = (dy_max / dy_min) - 1.0;
	az2 = (dz_max / dz_min) - 1.0;

	i_min = i = 0;
        x_grid[0] = x1;
        dx[0] = dx_min * stretch_factor(x_grid[0] - x2,ax1,ax2);
        while(1) {
                x_grid[i+1] = x_grid[i] + dx[i];
                if(x_grid[i+1] < x2) {
                        dx[i+1] = dx_min * stretch_factor(x_grid[i+1] - x2,ax1,ax2);
                }
                else if(x_grid[i+1] < x3) {
                        dx[i+1] = dx_min;
                }
                else if(x_grid[i+1] < x4) {
                        dx[i+1] = dx_min * stretch_factor(x_grid[i+1] - x3,ax1,ax2);
                }
                else if(x_grid[i+1] >= x4) {
                        i_max = i+1;
                        x_max = x_grid[i+1];
                        dx[i+1] = dx[i];        // not really needed
                        break;
                }
                ++i;
        }

	for(i = i_min; i <= i_max; i++) {
		fprintf(fp_dx,"%g\t%g\n",x_grid[i],dx[i]);
		fprintf(fp_dx_display,"%g\t%g\n",x_grid[i],y2);
		fprintf(fp_dx_display,"%g\t%g\n",x_grid[i],y3);
		fprintf(fp_x_grid_dat,"%d\t%g\t%g\n",i,x_grid[i],dx[i]);
	}

	j_min = j = 0;
	y_grid[0] = y_1;
	dy[0] = dy_min * stretch_factor(y_grid[0] - y2,ay1,ay2);
	while(1) {
		y_grid[j+1] = y_grid[j] + dy[j];
		if(y_grid[j+1] < y2) {
			dy[j+1] = dy_min * stretch_factor(y_grid[j+1] - y2,ay1,ay2);
		}
		else if(y_grid[j+1] < y3) {
			dy[j+1] = dy_min;
		}
		else if(y_grid[j+1] < y4) {
			dy[j+1] = dy_min * stretch_factor(y_grid[j+1] - y3,ay1,ay2);
		}
		else if(y_grid[j+1] >= y4) {
			j_max = j+1;
			y_max = y_grid[j+1];
			dy[j+1] = dy[j];	// not really needed
			break;
		}
		++j;
	}

	for(j = j_min; j <= j_max; j++) {
		fprintf(fp_dy,"%g\t%g\n",y_grid[j],dy[j]);
		fprintf(fp_dy_display,"%g\t%g\n",x2,y_grid[j]);
		fprintf(fp_dy_display,"%g\t%g\n",x3,y_grid[j]);
		fprintf(fp_y_grid_dat,"%d\t%g\t%g\n",j,y_grid[j],dy[j]);
	}

	k_min = k = 0;
	z_grid[0] = z1;
	dz[0] = dz_min * stretch_factor(z_grid[0] - z2,az1,az2);
	while(1) {
		z_grid[k+1] = z_grid[k] + dz[k];
		if(z_grid[k+1] < z2) {
			dz[k+1] = dz_min * stretch_factor(z_grid[k+1] - z2,az1,az2);
		}
		else if(z_grid[k+1] < z3) {
			dz[k+1] = dz_min;
		}
		else if(z_grid[k+1] < z4) {
			dz[k+1] = dz_min * stretch_factor(z_grid[k+1] - z3,az1,az2);
		}
		else if(z_grid[k+1] >= z4) {
			k_max = k+1;
			z_max = z_grid[k+1];
			dz[k+1] = dz[k];	// not really needed
			break;
		}
		++k;
	}

	for(k = k_min; k <= k_max; k++) {
		fprintf(fp_dz,"%g\t%g\n",z_grid[k],dz[k]);
		fprintf(fp_dz_display,"%g\t%g\n",x2,z_grid[k]);
		fprintf(fp_dz_display,"%g\t%g\n",x3,z_grid[k]);
		fprintf(fp_z_grid_dat,"%d\t%g\t%g\n",k,z_grid[k],dz[k]);
	}

	fprintf(stderr,"number of grid points = %d x %d x %d = %d\n",i_max+1,j_max+1,k_max+1,(i_max+1)*(j_max+1)*(k_max+1));
	fprintf(stderr,"grid points for uniform grid = %d x %d x %d = %d\n",(int)((x4-x1)/dx_min)+1,(int)((y4-y_1)/dy_min)+1,(int)((z4-z1)/dz_min)+1,((int)((x4-x1)/dx_min)+1)*((int)((y4-y_1)/dy_min)+1)*((int)((z4-z1)/dz_min)+1));
	fprintf(stderr,"x_max = %g\ty_max = %g\tz_max = %g\n",x_grid[i_max],y_grid[j_max],z_grid[k_max]);

	fclose(fp_dx);
	fclose(fp_dy);
	fclose(fp_dz);
	fclose(fp_dx_display);
	fclose(fp_dy_display);
	fclose(fp_dz_display);
	fclose(fp_x_grid_dat);
	fclose(fp_y_grid_dat);
	fclose(fp_z_grid_dat);

	return;
}

double
stretch_factor(double v, double a1, double a2)
{
	double val;
	double dum;
	double x;

	x = a1*fabs(v);
	dum = (exp(x) - exp(-x))/(exp(x) + exp(-x));
	val = 1.0 + a2 * pow(dum,2.0);

	return(val);
}		

	
place_recorder()
{
	FILE *fp_rec,*fp,*fp_u_init;
	int i,j,k,n,m,p,q;
        char tmp[N_Y_MAX],tmp2[N_Y_MAX];
        int i_rec,j_rec;
	int init_flag;
        char mes[800],c_junk;
	int n_layer[N_Z_MAX];
	int i_init_bare_min,i_init_bare_max;
	int j_init_bare_min,j_init_bare_max;
	int k_init_bare_min,k_init_bare_max;
	int im,jm,km;
	FILE *fp_marked;

	sprintf(mes,"%s/marked_points",output_directory);
	fp_marked = fopen(mes,"w");

        fp = fopen(recorder_geometry,"r");
        if(fp == NULL) {
                fprintf(stderr,"can't open %s for geometry data\n",recorder_geometry);
                exit(0);
        }

        while(1) {
                if(fgets(tmp,300,fp) == NULL) {
                        fprintf(stderr, "error 0 in reading geometry file\n");
                        (void)exit(0);
                }
                if(tmp[0] == '#') break;
        }

        fgets(tmp,300,fp);
        sscanf(tmp,"%lf",&dx_bare);
        fgets(tmp,300,fp);
        sscanf(tmp,"%d",&i_max_bare);
        fgets(tmp,300,fp);
        sscanf(tmp,"%d",&j_max_bare);
        fgets(tmp,300,fp);
        sscanf(tmp,"%d",&k_max_bare);
        fgets(tmp,300,fp);
        sscanf(tmp,"%d",&i_recorder_bare);
        fgets(tmp,300,fp);
        sscanf(tmp,"%d",&j_recorder_bare);
        fgets(tmp,300,fp);
        sscanf(tmp,"%d",&k_recorder_bare);

	x_recorder_bare = i_recorder_bare * dx_bare;
	y_recorder_bare = j_recorder_bare * dx_bare;
	z_recorder_bare = k_recorder_bare * dx_bare;

        for(i = 0; i <= i_max_bare; i++) {
                for(j = 0; j <= j_max_bare; j++) {
                	for(k = 0; k <= k_max_bare; k++) {
                        	bare_geometry[i][j][k] = OPEN;
			}
                }
        }
        for(i = 0; i <= i_max_bare; i++) {
                for(j = 0; j <= j_max_bare; j++) {
                	bare_geometry[i][j][0] = LOWER_Z_BOUNDARY;
                	bare_geometry[i][j][k_max_bare] = UPPER_Z_BOUNDARY;
		}
        }
        for(j = 0; j <= j_max_bare; j++) {
                for(k = 0; k <= k_max_bare; k++) {
                	bare_geometry[0][j][k] = LOWER_X_BOUNDARY;
                	bare_geometry[i_max_bare][j][k] = UPPER_X_BOUNDARY;
		}
        }
        for(i = 0; i <= i_max_bare; i++) {
                for(k = 0; k <= k_max_bare; k++) {
                	bare_geometry[i][0][k] = LOWER_Y_BOUNDARY;
                	bare_geometry[i][j_max_bare][k] = UPPER_Y_BOUNDARY;
		}
        }

        i = i_recorder_bare;
        j = j_recorder_bare;
        k = k_recorder_bare;	// k is the index of the next layer to be populated

// now read in the first "page" of the instrument from the rec_geom file
//	the pages are indexed using q starting with q = 1
// the velocity scale factor for each page will be stored in v_scale_factor_bare[q]

// each of these pages will be repeated a certain number of times creating
// the instrument layer by layer - each layer is parallel to the x-y plane

	q = 0;
        fgets(tmp,N_Y_MAX,fp);
	while(1) {
		if(tmp[0] != '@') {
            		fprintf(stderr,"error 1 in geometry file\n");
               		(void)exit(0);
        	}
		++q;
		sscanf(tmp+1,"%d %lf",&(n_layer[q]),&(v_scale_factor_bare[k]));	// n_layer[q] is the number of times to repeat this page
										// v_scale_factor_bare[k] is the scale factor for all of these layers
		if(n_layer[q] < 1) break;	// done with pages

fprintf(stderr,"page %d\tplaced starting in layer %d\trepeat %d times\tv_factor = %g\n",q,k,n_layer[q],v_scale_factor_bare[k]);
// starting a new page
// place the first layer of this page in the plane k 
        	i = i_recorder_bare;
        	while(1) {	// this loop reads in the page and places it in a single plane
               		if(fgets(tmp,N_Y_MAX,fp) == NULL) {
               	       		fprintf(stderr,"error 2 in geometry file\n");
               	        	(void)exit(0);
               	 	}
               	 	if(tmp[0] == '@') break;	// this means we are done with this page - start this loop over again for the next page

                	sscanf(tmp,"%d %s",&n,tmp2);	// n is the number of times to repeat this row (parallel to y direction)
							// tmp2 contains the data for this row
fprintf(stderr,"%s\n",tmp2);
fprintf(stderr,"now repeat this row %d times\n",n);
                	for(m = 0; m < n; m++) {
                       		j = j_recorder_bare;
                        	for(p = 0; p < strlen(tmp2); p++) {	// read the data for this row one value of j (y direction) at a time
                                	if(tmp2[p] == 'X') bare_geometry[i][j][k] = SOLID;
                                	if(tmp2[p] == '.') bare_geometry[i][j][k] = OPEN;
                                	if(tmp2[p] == 'I') {
						bare_geometry[i][j][k] = INIT_POINT;
//fprintf(stderr,"init point at %d\t%d\t%d\tv_scale = %g\n",i,j,k,v_scale_factor_bare[k]);
//fprintf(fp_diag,"%d\t%d\t%d\t%d\n",i,j,k,geometry[i][j][k]);
					}
                                	if(tmp2[p] == 'M') {	// these are special marking points
						bare_geometry[i][j][k] = MARKED;
					}
					++j;
				}
				++i;
			}
		}
// now repeat this layer n_layer[q]-1 times - the layer is in the x-y plane
	        for(p = 1; p < n_layer[q]; p++) {
			for(i = 0; i <= i_max_bare; i++) {
				for(j = 0; j <= j_max_bare; j++) {
					bare_geometry[i][j][k + p] = bare_geometry[i][j][k];
				}
			}
			v_scale_factor_bare[k + p] = v_scale_factor_bare[k];
		}  
		k += n_layer[q];
	}

// now find the initialization points in the channel
	i_init_bare_min = i_max_bare;
	i_init_bare_max = 0;
	j_init_bare_min = j_max_bare;
	j_init_bare_max = 0;
	k_init_bare_min = k_max_bare;
	k_init_bare_max = 0;
fprintf(fp_diag,"init points in bare and final coords\n");
fprintf(stderr,"init points in bare and final coords\n");
        for(i = 0; i <= i_max_bare; i++) {
                for(j = 0; j <= j_max_bare; j++) {
                	for(k = 0; k <= k_max_bare; k++) {
                        	if(bare_geometry[i][j][k] == INIT_POINT) {
					if(i < i_init_bare_min) i_init_bare_min = i;
					if(i > i_init_bare_max) i_init_bare_max = i;
					if(j < j_init_bare_min) j_init_bare_min = j;
					if(j > j_init_bare_max) j_init_bare_max = j;
					if(k < k_init_bare_min) k_init_bare_min = k;
					if(k > k_init_bare_max) k_init_bare_max = k;
//fprintf(fp_diag,"%d\t%d\t%d\t%d\t%d\t%d\tv_scale = %g\n",i,i_trans(i),j,j_trans(j),k,k_trans(k),v_scale_factor_bare[k]);

				}
			}
                }
        }
fprintf(fp_diag,"done with init points in bare coords\n");
fprintf(stderr,"done with init points in bare coords\n");

	i_init_min = i_trans(i_init_bare_min);
	i_init_max = i_trans(i_init_bare_max);
	j_init_min = j_trans(j_init_bare_min);
	j_init_max = j_trans(j_init_bare_max);
	k_init_min = k_trans(k_init_bare_min);
	k_init_max = k_trans(k_init_bare_max);

fprintf(fp_diag,"init points in channel\n");
fprintf(fp_diag,"%d-%d\t%d-%d\t%d-%d\n",
i_init_bare_min,i_init_bare_max,j_init_bare_min,j_init_bare_max,k_init_bare_min,k_init_bare_max);

// the next line contains info for u_init.dat
	fgets(tmp,300,fp);
	sprintf(tmp2,"%s/u_init.dat",output_directory);
	fp_u_init = fopen(tmp2,"w");
	fprintf(fp_u_init,"%s\n",tmp);
	fclose(fp_u_init);

// places where sound is recorded
/*
        fscanf(fp,"%d %d %d",&i_record_1_bare,&j_record_1_bare,&k_record_1_bare);
        fscanf(fp,"%d %d %d",&i_record_2_bare,&j_record_2_bare,&k_record_2_bare);
        fscanf(fp,"%d %d %d",&i_record_3_bare,&j_record_3_bare,&k_record_3_bare);
*/
	n_record = 1;
	while(1) {
		if(fscanf(fp,"%d %d %d",&(i_record_bare[n_record]),&(j_record_bare[n_record]),&(k_record_bare[n_record])) == EOF) break;
fprintf(stderr,"%d %d %d\n",i_record_bare[n_record],j_record_bare[n_record],k_record_bare[n_record]);
		++n_record;
		if(n_record > 10) break;
	}
	--n_record;

        fclose(fp);

//	now convert from the bare geometry (uniform grid) to the new geometry (nonuniform grid)
//	will have to adjust the placement of some points but keep them close to the original
//	locations in terms of real space (x and y)

//	first locate the grid point closest to the desired location of the lower left
//	point on the recorder

	for(i = 0; i < i_max; i++) {
		if(fabs(x_grid[i] - x_recorder_bare) <= dx[i]/2.0) {
			i_recorder_corner = i;
			break;
		}
	}
	for(j = 0; j < j_max; j++) {
		if(fabs(y_grid[j] - y_recorder_bare) <= dy[j]/2.0) {
			j_recorder_corner = j;
			break;
		}
	}
	for(k = 0; k < k_max; k++) {
		if(fabs(z_grid[k] - z_recorder_bare) <= dz[k]/2.0) {
			k_recorder_corner = k;
			break;
		}
	}

//	now transfer the recorder location in bare_geometry[][][] to points in geometry[][][]

	for(i = 0; i < i_max; i++) {
		for(j = 0; j < j_max; j++) {
			for(k = 0; k < k_max; k++) {
				geometry[i][j][k] = OPEN;
			}
		}
	}
	for(i = 0; i <= i_max; i++) {
		for(k = 0; k < k_max; k++) {
                	geometry[i][0][k] = LOWER_Y_BOUNDARY;
                	geometry[i][j_max][k] = UPPER_Y_BOUNDARY;
		}
        }
        for(j = 0; j <= j_max; j++) {
		for(k = 0; k < k_max; k++) {
                	geometry[0][j][k] = LOWER_X_BOUNDARY;
                	geometry[i_max][j][k] = UPPER_X_BOUNDARY;
		}
        }
        for(i = 0; i <= i_max; i++) {
        	for(j = 0; j <= j_max; j++) {
                		geometry[i][j][0] = LOWER_Z_BOUNDARY;
                		geometry[i][j][k_max] = UPPER_Z_BOUNDARY;
		}
        }

fprintf(fp_diag,"converting\n");
fprintf(stderr,"converting\n");
/*	this is the old version - note that the inner-most loop can overwrite points that might have been INIT_POINTS in previous loops 
	for(i = 0; i < i_max_bare; i++) {
		for(j = 0; j < j_max_bare; j++) {
			for(k = 0; k < k_max_bare; k++) {
				if(bare_geometry[i][j][k] == SOLID) {
					geometry[i_trans(i)][j_trans(j)][k_trans(k)] = SOLID;
				}
				else if(bare_geometry[i][j][k] == INIT_POINT) {
					geometry[i_trans(i)][j_trans(j)][k_trans(k)] = INIT_POINT;
fprintf(fp_diag,"%d\t%d\t%d\t%d\t%d\t%d\n",i,i_trans(i),j,j_trans(j),k,k_trans(k));
				}
			}
		}
	}
*/
// this is the new version - this should get rid of over-writing problem noted above
// deal with marked points first
	for(i = 0; i < i_max_bare; i++) {
		for(j = 0; j < j_max_bare; j++) {
			for(k = 0; k < k_max_bare; k++) {
				if(bare_geometry[i][j][k] == MARKED) {
					bare_geometry[i][j][k] = OPEN;
					im = i_trans(i);
					jm = j_trans(j);
					km = k_trans(k);
					fprintf(fp_marked,"%d\t%d\t%d\t%g\t%g\t%g\n",im,jm,km,x_grid[im],y_grid[jm],z_grid[km]);
// fprintf(stderr,"found a marked point\n");
// exit(0);
				}
			}
		}
	}
	fclose(fp_marked);

	for(i = 0; i < i_max_bare; i++) {
		for(j = 0; j < j_max_bare; j++) {
			for(k = 0; k < k_max_bare; k++) {
				if(bare_geometry[i][j][k] == SOLID) {
					geometry[i_trans(i)][j_trans(j)][k_trans(k)] = SOLID;
				}
			}
		}
	}
	for(i = 0; i < i_max_bare; i++) {
		for(j = 0; j < j_max_bare; j++) {
			for(k = 0; k < k_max_bare; k++) {
				if(bare_geometry[i][j][k] == INIT_POINT) {
					geometry[i_trans(i)][j_trans(j)][k_trans(k)] = INIT_POINT;
					v_scale_factor[k_trans(k)] = v_scale_factor_bare[k];
fprintf(fp_diag,"%d\t%d\t%d\t%d\t%d\t%d\n",i,i_trans(i),j,j_trans(j),k,k_trans(k));
				}
			}
		}
	}
fprintf(fp_diag,"finished converting\n");
fprintf(stderr,"finished converting\n");

/*
	i_init_min = i_trans(i_init_min_bare);
	i_init_max = i_trans(i_init_max_bare);
	j_init_min = j_trans(j_init_min_bare);
	j_init_max = j_trans(j_init_max_bare);
	k_init_min = k_trans(k_init_min_bare);
	k_init_max = k_trans(k_init_max_bare);
*/

fprintf(fp_diag,"checking\n");
for(i = i_init_min; i <= i_init_max; i++) {
	for(j = j_init_min; j <= j_init_max; j++) {
		for(k = k_init_min; k <= k_init_max; k++) {
			if(geometry[i][j][k] != INIT_POINT) {
fprintf(fp_diag,"%d\t%d\t%d\n",i,j,k);
			}
		}
	}
}
fprintf(fp_diag,"done checking\n");
fprintf(stderr,"done checking\n");
fprintf(stderr,"n_record = %d\n",n_record);

	for(i = 1; i <= n_record; i++) {
		i_record[i] = i_trans(i_record_bare[i]);
		j_record[i] = j_trans(j_record_bare[i]);
		k_record[i] = k_trans(k_record_bare[i]);
fprintf(stderr,"%d\ti/j/k:%d / %d / %d\n",i,i_record[i],j_record[i],k_record[i]);
	}

/*
        i_record_1 = i_trans(i_record_1_bare);
	j_record_1 = j_trans(j_record_1_bare);
	k_record_1 = k_trans(k_record_1_bare);
        i_record_2 = i_trans(i_record_2_bare);
	j_record_2 = j_trans(j_record_2_bare);
	k_record_2 = k_trans(k_record_2_bare);
        i_record_3 = i_trans(i_record_3_bare);
	j_record_3 = j_trans(j_record_3_bare);
	k_record_3 = k_trans(k_record_3_bare);
*/
        
        return;
}

/*	
 *	functions to convert from bare coordinates to coordinates on the nonuniform grid
 */

i_trans(i_bare)
int i_bare;
{
	int i;
	double x_bare;
	static int i_last = 0;

	x_bare = i_bare * dx_bare;

	if(fabs(x_bare - x_grid[0]) <= dx[0]/2.0) {
		i_last = 0;
		return(i_last);
	}
	for(i = i_last; i < i_max; i++) {
		if(fabs(x_grid[i] - x_bare) <= dx[i-1]/2.0) {
			i_last = i;
			return(i_last);
		}
		else if(fabs(x_bare - x_grid[i]) <= dx[i]/2.0) {
			i_last = i;
			return(i_last);
		}
	}
	if(fabs(x_grid[i_max] - x_bare) <= dx[i_max-1]/2.0) {
		i_last = i_max;
		return(i_last);
	}
	if(fabs(x_bare - x_grid[i_max]) <= dx[i_max]/2.0) {
		i_last = i_max;
		return(i_last);
	}
	for(i = 0; i < i_last+1; i++) {
		if(fabs(x_grid[i] - x_bare) <= dx[i-1]/2.0) {
			i_last = i;
			return(i_last);
		}
		else if(fabs(x_bare - x_grid[i]) <= dx[i]/2.0) {
			i_last = i;
			return(i_last);
		}
	}
	fprintf(stderr,"can't locate the new x coordinate\n");
	fprintf(stderr,"i_bare = %d\tx_bare = %g\n",i_bare,x_bare);
	fprintf(fp_log,"can't locate the new x coordinate\n");
	fprintf(fp_log,"i_bare = %d\tx_bare = %g\n",i_bare,x_bare);
	fflush(fp_log);
	exit(0);
	
	return(i_last);
}

j_trans(j_bare)
int j_bare;
{
	int j;
	double y_bare;
	static int j_last = 0;

	y_bare = j_bare * dx_bare;	// bare coordinate of the corner of the location to be converted

	if(fabs(y_bare - y_grid[0]) <= dy[0]/2.0) {
		j_last = 0;
		return(j_last);
	}
	for(j = j_last; j < j_max; j++) {
		if(fabs(y_grid[j] - y_bare) <= dy[j-1]/2.0) {
			j_last = j;
			return(j_last);
		}
		else if(fabs(y_bare - y_grid[j]) <= dy[j]/2.0) {
			j_last = j;
			return(j_last);
		}
	}
	if(fabs(y_grid[j_max] - y_bare) <= dy[j_max-1]/2.0) {
		j_last = j_max;
		return(j_last);
	}
	if(fabs(y_bare - y_grid[j_max]) <= dy[j_max]/2.0) {
		j_last = j_max;
		return(j_last);
	}
	for(j = 1; j < j_last+1; j++) {
		if(fabs(y_grid[j] - y_bare) <= dy[j-1]/2.0) {
			j_last = j;
			return(j_last);
		}
		else if(fabs(y_bare - y_grid[j]) <= dy[j]/2.0) {
			j_last = j;
			return(j_last);
		}
	}
	fprintf(stderr,"can't locate the new y coordinate\n");
	fprintf(stderr,"j_bare = %d\ty_bare = %g\n",j_bare,y_bare);
	fprintf(fp_log,"can't locate the new y coordinate\n");
	fprintf(fp_log,"j_bare = %d\ty_bare = %g\n",j_bare,y_bare);
	fflush(fp_log);
	exit(0);
	
	return(0);
}

k_trans(k_bare)
int k_bare;
{
	int k;
	double z_bare;
	static int k_last = 0;

	z_bare = k_bare * dx_bare;

	if(fabs(z_bare - z_grid[0]) <= dz[0]/2.0) {
		k_last = 0;
		return(k_last);
	}
	for(k = k_last; k < k_max; k++) {
		if(fabs(z_grid[k] - z_bare) <= dz[k-1]/2.0) {
			k_last = k;
			return(k_last);
		}
		else if(fabs(z_bare - z_grid[k]) <= dz[k]/2.0) {
			k_last = k;
			return(k_last);
		}
	}
	if(fabs(z_grid[k_max] - z_bare) <= dz[k_max-1]/2.0) {
		k_last = k_max;
		return(k_last);
	}
	if(fabs(z_bare - z_grid[k_max]) <= dz[k_max]/2.0) {
		k_last = k_max;
		return(k_last);
	}
	for(k = 1; k <= k_last; k++) {
		if(fabs(z_grid[k] - z_bare) <= dz[k-1]/2.0) {
			k_last = k;
			return(k_last);
		}
		else if(fabs(z_bare - z_grid[k]) <= dz[k]/2.0) {
			k_last = k;
			return(k_last);
		}
	}
	fprintf(stderr,"can't locate the new x coordinate\n");
	fprintf(stderr,"k_bare = %d\tz_bare = %g\n",k_bare,z_bare);
	fprintf(fp_log,"can't locate the new z coordinate\n");
	fprintf(fp_log,"k_bare = %d\tz_bare = %g\n",k_bare,z_bare);
	fflush(fp_log);
	exit(0);
	
	return(0);
}

/*
 *	store geometry data in a vtk file for view with paraview
 */
make_vtk_file1()
{
	int i,j,k;
	FILE *fp_vtk;
	long int n_points;
	int n_rec_points;
	int i_rec_min,j_rec_min,k_rec_min;
	int i_rec_max,j_rec_max,k_rec_max;
	char mes[200];

	i_rec_min = i_max;
	j_rec_min = j_max;
	k_rec_min = k_max;
	i_rec_max = j_rec_max = k_rec_max = 0;
	for(i = 0; i <= i_max; i++) {
		for(j = 0; j <= j_max; j++) {
			for(k = 0; k <= k_max; k++) {
				if((geometry[i][j][k] == SOLID) 
					|| (geometry[i][j][k] == INIT_POINT)) {
					if(i < i_rec_min) i_rec_min = i;
					if(i > i_rec_max) i_rec_max = i;
					if(j < j_rec_min) j_rec_min = j;
					if(j > j_rec_max) j_rec_max = j;
					if(k < k_rec_min) k_rec_min = k;
					if(k > k_rec_max) k_rec_max = k;
				}
			}
		}
	}

	n_points = (i_rec_max-i_rec_min+1)*(j_rec_max-j_rec_min+1)*(k_rec_max-k_rec_min+1);

	sprintf(mes,"%s/geometry.vtk",vtk_directory);
	fp_vtk = fopen(mes,"w");

	fprintf(fp_vtk,"# vtk DataFile Version 2.0\n");
	fprintf(fp_vtk,"Recorder Dataset\n");
	fprintf(fp_vtk,"ASCII\n");
	fprintf(fp_vtk,"DATASET STRUCTURED_GRID\n");
	fprintf(fp_vtk,"DIMENSIONS %d %d %d\n",k_rec_max-k_rec_min+1,j_rec_max-j_rec_min+1
									,i_rec_max-i_rec_min+1);
									
	fprintf(fp_vtk,"POINTS %ld float\n",n_points);

	for(i = i_rec_min; i <= i_rec_max; i++) {
		for(j = j_rec_min; j <= j_rec_max; j++) {
			for(k = k_rec_min; k <= k_rec_max; k++) {
//				fprintf(fp_vtk,"%g\t%g\t%g\n",dx_min*i,dy_min*j,dz_min*k);
				fprintf(fp_vtk,"%g\t%g\t%g\n",x_grid[i],y_grid[j],z_grid[k]);
			}
		}
	}

	fprintf(fp_vtk,"\n");
	fprintf(fp_vtk,"POINT_DATA %ld\n",n_points);
	fprintf(fp_vtk,"SCALARS rho float 1\n");
	fprintf(fp_vtk,"LOOKUP_TABLE default\n");

	for(i = i_rec_min; i <= i_rec_max; i++) {
		for(j = j_rec_min; j <= j_rec_max; j++) {
			for(k = k_rec_min; k <= k_rec_max; k++) {
				if(geometry[i][j][k] == SOLID) {
					fprintf(fp_vtk,"10.0\n");
				}
				else if(geometry[i][j][k] == INIT_POINT) {
					fprintf(fp_vtk,"5.0\n");
				}
//				else if((i == i_record_1) && (j == j_record_1) && (k == k_record_1)) {
//					fprintf(fp_vtk,"0.25\n");
//				}
//				else if((i == i_record_2) && (j == j_record_2) && (k == k_record_2)) {
//					fprintf(fp_vtk,"0.25\n");
//				}
//				else if((i == i_record_3) && (j == j_record_3) && (k == k_record_3)) {
//					fprintf(fp_vtk,"0.25\n");
//				}
				else {
					fprintf(fp_vtk,"0.0\n");
				}
			}
		}
	}

	fclose(fp_vtk);

	return;
}


/*
 *      save the geometry of the simulation:
 *	the solid portion of the recorder
 *      the points where the velocity is initialized (and held fixed) 
 *      the points where rho-vs-time is recordeded 
 */
save_geometry()
{
        int i,j,k;
        FILE *fp_solid,*fp_initialize,*fp_sound,*fp_instrument;
        char mes[400],tmp[200];

fprintf(stderr,"starting to save geom\n");

	sprintf(mes,"%s/geometry_solid.dat",output_directory);
        fp_solid = fopen(mes,"w");
	sprintf(mes,"%s/geometry_initialize.dat",output_directory);
        fp_initialize = fopen(mes,"w");
	sprintf(mes,"%s/geometry_sound.dat",output_directory);
        fp_sound = fopen(mes,"w");
	sprintf(mes,"%s/geometry_instrument.dat",output_directory);
        fp_instrument = fopen(mes,"w");

fprintf(fp_diag,"saving geom\n");
fprintf(stderr,"saving geom\n");
        for(i = 0; i <= i_max; i++) {
                for(j = 0; j <= j_max; j++) {
                	for(k = 0; k <= k_max; k++) {
                       		if(geometry[i][j][k] == INIT_POINT) {
                               		fprintf(fp_initialize,"%d\t%d\t%d\t%g\n",i,j,k,v_scale_factor[k]);
fprintf(fp_diag,"%d\t%d\t%d\t%d\t%g\n",i,j,k,geometry[i][j][k],v_scale_factor[k]);
				}
                        	else if(geometry[i][j][k] !=  OPEN) {
                                	fprintf(fp_solid,"%d\t%d\t%d\n",i,j,k);
				}
                        }
		}
	}
        for(i = 1; i < i_max; i++) {
                for(j = 1; j < j_max; j++) {
                	for(k = 1; k < k_max; k++) {
                        	if(geometry[i][j][k] !=  OPEN) {
                                	fprintf(fp_instrument,"%d\t%d\t%d\n",i,j,k);
				}
                        }
		}
	}

fprintf(fp_diag,"finished saving geom\n");
fflush(fp_diag);
fprintf(stderr,"finished saving geom\n");
fprintf(stderr,"n_record = %d\n",n_record);
	for(i = 1; i <= n_record; i++) {
		fprintf(fp_sound,"%d\t%d\t%d\n",i_record[i],j_record[i],k_record[i]);
fprintf(stderr,"%d\t%d\t%d\n",i_record[i],j_record[i],k_record[i]);
	}

/*
	fprintf(fp_sound,"%d\t%d\t%d\n",i_record_1,j_record_1,k_record_1);
	fprintf(fp_sound,"%d\t%d\t%d\n",i_record_2,j_record_2,k_record_2);
	fprintf(fp_sound,"%d\t%d\t%d\n",i_record_3,j_record_3,k_record_3);
*/

	fclose(fp_solid);
	fclose(fp_initialize);
	fclose(fp_sound);
	fclose(fp_instrument);

	return;
}

/*
 *	store bare geometry data in a vtk file for view with paraview
 */
make_bare_vtk_file()
{
	int i,j,k;
	FILE *fp_vtk;
	long int n_points;
	int n_rec_points;
	int i_rec_min,j_rec_min,k_rec_min;
	int i_rec_max,j_rec_max,k_rec_max;
	char mes[200];

	i_rec_min = i_max_bare;
	j_rec_min = j_max_bare;
	k_rec_min = k_max_bare;
	i_rec_max = j_rec_max = k_rec_max = 0;
	for(i = 0; i <= i_max; i++) {
		for(j = 0; j <= j_max; j++) {
			for(k = 0; k <= k_max; k++) {
				if((bare_geometry[i][j][k] == SOLID) 
					|| (bare_geometry[i][j][k] == INIT_POINT)) {
					if(i < i_rec_min) i_rec_min = i;
					if(i > i_rec_max) i_rec_max = i;
					if(j < j_rec_min) j_rec_min = j;
					if(j > j_rec_max) j_rec_max = j;
					if(k < k_rec_min) k_rec_min = k;
					if(k > k_rec_max) k_rec_max = k;
				}
			}
		}
	}

	n_points = (i_rec_max-i_rec_min+1)*(j_rec_max-j_rec_min+1)*(k_rec_max-k_rec_min+1);

	sprintf(mes,"%s/bare_geometry.vtk",vtk_directory);
	fp_vtk = fopen(mes,"w");

	fprintf(fp_vtk,"# vtk DataFile Version 2.0\n");
	fprintf(fp_vtk,"Recorder Dataset\n");
	fprintf(fp_vtk,"ASCII\n");
	fprintf(fp_vtk,"DATASET STRUCTURED_GRID\n");
	fprintf(fp_vtk,"DIMENSIONS %d %d %d\n",k_rec_max-k_rec_min+1,j_rec_max-j_rec_min+1
									,i_rec_max-i_rec_min+1);
									
	fprintf(fp_vtk,"POINTS %ld float\n",n_points);

	for(i = i_rec_min; i <= i_rec_max; i++) {
		for(j = j_rec_min; j <= j_rec_max; j++) {
			for(k = k_rec_min; k <= k_rec_max; k++) {
				fprintf(fp_vtk,"%g\t%g\t%g\n",dx_min*i,dy_min*j,dz_min*k);
			}
		}
	}

	fprintf(fp_vtk,"\n");
	fprintf(fp_vtk,"POINT_DATA %ld\n",n_points);
	fprintf(fp_vtk,"SCALARS rho float 1\n");
	fprintf(fp_vtk,"LOOKUP_TABLE default\n");

	for(i = i_rec_min; i <= i_rec_max; i++) {
		for(j = j_rec_min; j <= j_rec_max; j++) {
			for(k = k_rec_min; k <= k_rec_max; k++) {
				if(bare_geometry[i][j][k] == SOLID) {
					fprintf(fp_vtk,"10.0\n");
				}
				else if(bare_geometry[i][j][k] == INIT_POINT) {
					fprintf(fp_vtk,"5.0\n");
				}
				else {
					fprintf(fp_vtk,"0.0\n");
				}
			}
		}
	}

	fclose(fp_vtk);

	return;
}

check_init_region()
{
	int i,j,k;

	fprintf(fp_diag,"checking init region\n");
	for(i = i_init_min; i <= i_init_max; i++) {
		for(j = j_init_min; j <= j_init_max; j++) {
			for(k = k_init_min; k <= k_init_max; k++) {
//				if(geometry[i][j][k] != INIT_POINT) {
					fprintf(fp_diag,"%d\t%d\t%d\t%d\n",i,j,k,geometry[i][j][k]);
//				}
			}
		}
	}
	fprintf(fp_diag,"done checking init region\n");

	return;
}

/*
 *	store geometry data in a second vtk file for viewing with paraview
 *	this file has only the solid points
 */
make_vtk_file2()
{
	int i,j,k;
	FILE *fp_vtk;
	long int n_points;
	int n_rec_points;
	int i_rec_min,j_rec_min,k_rec_min;
	int i_rec_max,j_rec_max,k_rec_max;
	char mes[200];

	i_rec_min = i_max;
	j_rec_min = j_max;
	k_rec_min = k_max;
	i_rec_max = j_rec_max = k_rec_max = 0;
	n_points = 0;
	for(i = 0; i <= i_max; i++) {
		for(j = 0; j <= j_max; j++) {
			for(k = 0; k <= k_max; k++) {
				if(geometry[i][j][k] == SOLID) {
					++n_points; 
				}
			}
		}
	}

	sprintf(mes,"%s/geometry_2.vtk",vtk_directory);
	fp_vtk = fopen(mes,"w");

	fprintf(fp_vtk,"# vtk DataFile Version 2.0\n");
	fprintf(fp_vtk,"Recorder Dataset\n");
	fprintf(fp_vtk,"ASCII\n");
	fprintf(fp_vtk,"DATASET UNSTRUCTURED_GRID\n");
//	fprintf(fp_vtk,"DIMENSIONS %d %d %d\n",k_rec_max-k_rec_min+1,j_rec_max-j_rec_min+1,i_rec_max-i_rec_min+1);
									
	fprintf(fp_vtk,"POINTS %ld float\n",n_points);

	for(i = 0; i <= i_max; i++) {
		for(j = 0; j <= j_max; j++) {
			for(k = 0; k <= k_max; k++) {
				if(geometry[i][j][k] == SOLID) 	
					 fprintf(fp_vtk,"%g\t%g\t%g\n",x_grid[i],y_grid[j],z_grid[k]);
			}
		}
	}

	fprintf(fp_vtk,"\n");
	fprintf(fp_vtk,"POINT_DATA %ld\n",n_points);
	fprintf(fp_vtk,"SCALARS rho float 1\n");
	fprintf(fp_vtk,"LOOKUP_TABLE default\n");

	for(i = 0; i <= i_max; i++) {
		for(j = 0; j <= j_max; j++) {
			for(k = 0; k <= k_max; k++) {
				if(geometry[i][j][k] == SOLID) {
					fprintf(fp_vtk,"1.0\n");
				}
			}
		}
	}

	fclose(fp_vtk);

	return;
}
